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Abstract 

We propose a new constrained level-set method for semi-automatic image 
segmentation. The method allows to specify which parts of the image lie 
inside respectively outside the segmented objects. Such an a-priori infor- 
mation can be expressed in terms of upper and lower constraints prescribed 
for the level-set function. Constraints have the same meaning as the initial 
seeds of the graph-cuts based methods for image segmentation. A numerical 
approximation scheme is based on the complementary-finite volumes method 
combined with the Projected successive over-relaxation method adopted for 
solving range bounds constrained linear complementarity problems. The 
advantage of the constrained level-set method is demonstrated on several 
artificial images as well as on cardiac MRI data. 

Keywords: linear complementarity, image processing, segmentation, 
level-set method, projected successive over-relaxation method, graph cuts 
2000 MSC: 90C33, 35K55, 35K52, 53C44, 74S10, 74G15 



1. Introduction 

The level-set methods for the image segmentation have been studied and 
applied during the last two decades. The level-set method applied in the im- 
age segmentation is typically an iterative method. The segmentation starts 
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with an initial curve representing an initial guess for the segmented ob- 
ject and it is evolved towards the segmented object by means of a suitable 
geometric law taking into account the direction to the segmented object and 
also the curvature of evolved curves. Loosely speaking, the better the initial 
guess is, the better and faster the segmentation process is. This is profitable 
for processing of time sequences where the final segmentation of one frame 
may serve as the initial guess for the next frame. We refer the reader to a 
wide range of literature on this topic like Caselles et al. [i5|], Handlovicova 



et al. j8||, Osher, Paragios [14] or Sethian ll5|] and references therein. In 
comparison to parametric models studied e.g. in Benes et al. [l| and Kass 
et al. [9] the level-set methods can handle topological changes and so one 
initial curve can split and segment more separate objects. 

Very different segmentation methods are the graph-cuts methods (see. 
Boykov et al. (2I, ^) which are based on graph theory and algorithms for 
finding the minimal cuts respectively the maximal flows. It is also an op- 
timization problem but algorithms for the minimal-cuts/maximal-fiows find 
globally the best solution which is usually not the case of the level-set meth- 
ods. These algorithms are not iterative and they do not require initial curves. 
Instead of it, they need initial seeds - one or more points or lines in the interior 
and exterior of the segmented object. 

In this article, we show how to incorporate the a-priori information in 
form of the initial seeds to the level-set method. Our main purpose is to 
propose a new constrained level-set method which can be applied to the image 
segmentation problems. In comparison to the classical level-set methods 
(c.f. 0), our method allows an expert to prescribe an a-priori information 
by marking parts which are surely inside or outside the segmented region. 
Such a possibility can be helpful for an expert conducting an interactive 
segmentation of medical data. 

A numerical approximation scheme is based on the complementary-finite 
volumes method developed by Handlovicova et al. in fsl] combined with the 
Projected successive over-relaxation method for solving constrained problems 
proposed by Mangasarian in and Elliott, Ockendon in [7]. The advan- 
tage of the constrained level-set method is demonstrated on several artificial 
images as well as on cardiac MRI data. 

The paper is organized as follows. In Section [2] we show the common 
level-set method for the image segmentation together with the numerical 
scheme and successive over-relaxation (SOR) method. Section [3] explains 
the constrained level-set method with appropriate numerical scheme. As a 
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solver for the linear complementarity problem with range bounds we adopt 
the Projected SOR method. Comparison with the common level-set method 
and contributions of the constrained level-set method are demonstrated in 
Section [H 



2. The level-set method for the image segmentation 

2.1. Time-space continuous framework of the level-set method 

We assume an image is represented by the greyscale image function Iq : 
Q [0, 1] defined on a two dimensional rectangle Q = [0, Li] x [0,L2]. A 
common idea how to segment an object in the image is to start from a 
closed, embedded and smooth initial curve approximating the shape of 
the object and let it evolve towards the exact boundary of the object. To 
this end, we construct a family of evolving curves with the property that 
converges to the boundary of a segmented object as t goes to infinity. 
There are many ways how to construct such a fiow of planar curves. Among 
them we will focus our attention to the fiow of curves proposed in the active 
contour model (c.f. Caselles et al. in [6], Kichenassamy et al. in [10]). A 
problem of finding a boundary of an object in the image can be reformulated 
as a problem of construction of planar curves on which the gradient V/q of 
the image intensity function Iq is large. 

Assuming is a smooth curve we can evaluate the unit tangent vector 
T(x) and outer unit normal vector N(x). Each point x G ^Ms evolved in the 
normal direction with the normal velocity l/(x, t) = 5tX • N(x, t). Although 
the velocity vector ^^x can be decomposed into its tangential and normal 
parts, it should be noted that only the motion in the normal direction has 
impact on the shape of the closed curve G^- 

Following the active contour model (c.f. jsl, lo|), Mikula and Sevcovic 



13| considered a generalized form of the normal velocity: 

V (x, t) = (x) H (x, t) + Vg' (x) ■ N (x, t) , (1) 

where H is the curvature of and = g {\Ga V/o|). Here g is Sl smooth 
edge detector function g : [0, oc) (0, oc) such that g' < 0, ^'(0) = 1, g{+oo) = 
and g\s) < Cg{s)^ < C, 5 > 0, for some constant (7 > 0. A typical 

example is the function g{s) = 1/(1 + Xs'^) where A > is the contrast param- 
eter. Notice that, for a given smooth intensity function /q, the vector field 
W{x) = —Vg^{x) has an important geometric property as it points towards 
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edges in the image where the norm of the gradient V/q is large (c.f. fl3^). 
Notice that a possible lack of smoothness of Iq (e.g. due to a noise) can be 
overcome by taking the convolution of Iq with a smooth Gaussian moUifier 
with the dispersion > (see (lil, fl2j). The term (x) • ^ h^t) 
pushes the evolved curve Q towards the edge of the image Iq 
The effect of the curvature term H (x, t) consists in smoothing the segmented 
curve by means of minimization of its total length. This property makes the 
segmentation model robust for application even in the case of a noisy image. 
Notice that the term slows down the normal velocity in the vicinity of an 
edges of Iq (c.f. In the level-set method is given implicitly as 

= {^en\u{^,t) = 0} , 

where 2/ is a real valued smooth function defined on Q such that u{x.) < 
for all X belonging to the interior of and (x) > for all x belonging 
to the exterior of a Jordan curve G^- Following derivation from [12^, the 
level-set formulation of ([T]) can be stated in terms of a solution u to the 
initial-boundary value problem 

Ut = QV-^g'^) mQx{0,Ti (2) 
d^u = on dn, (3) 

u\t=o= Uini in f], (4) 

where Uini is the initial level-set function corresponding to the initial curve 
^ d^u = Vu • u and u is the outer normal unit vector of the boundary 
dQ of a computational domain The quantity Q should be equal to \\/u\. 
However, for practical purposes, it is regularized by means of the Tichonov 
regularization, i.e. 



Q = ^Je^ + \Vu\\ (5) 
where < e <C 1 is a small regularizing parameter. 

2.2. Time-space discretized framework of the level-set method 

We discretize the initial-boundary value problem ([2j)-(|lj) by means of 
the method of complementary finite volumes developed by Handlovicova et 
al. Q in the context of a class of level-set equations arising in the image 
processing. Let r be a time step for time discretization. Let h = (/ii, /12) be 
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spatial discretization steps such that hi = jf- for some A^^ G N+, i = 1,2. We 
define a numerical grid, its closure and its boundary as 

c^h = {{ihujh2)\i = l"'Ni-l,j = 1---N2-1}, 
uJh = {{ihujh2)\i = 0"'NiJ = 0---N2}, 



For a function u G C {Q x {O^T]] R) we define its piecewise constant approxi- 
mation on uJ at the time kr as a grid function defined by u^j = u (i/ii, J'7i2, kr). 
We furthermore introduce a dual mesh Vh defined as 



X 



l...iVi-l,j = l...iV2-l 



For < i < A^i, < j < A^2, ^ and j fixed, we consider a finite volume Vij 
of the dual mesh Vh- We denote its interior by its boundary by Vij and 
let lJi{vij) be the volume of Vtij (see Fig. [1]). We also denote the set of all 
neighboring volumes of a finite volume Vij by Mij. For all inner finite volumes 
Vij of the dual mesh V/^, the boundary Vij consists of four linear segments. 
We denote them as F^^ - for ij G Nij. It means that F^^ - is a boundary 
between the finite volumes Vij and v^j. Let /^^ - be the length of this part of 
Vij. Dividing equation ([2j) by the term integrating it over the interior 
of a finite volume Vij we obtain the equation: 

Applying the Euler backward in time discretization of Ut we end up with the 
following time semi-discretization of (Ej): 



Next, after long and tedious but straightforward calculations and evaluations 
described in Appendix we are in a position to formulate full time-space dis- 
cretization of the level set equation ([2j) . It can be rewritten in a form of the 
following system of linear equations: 
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Figure 1: A description of a finite volume Vij defined on the dual mesh. Here Qij is its 
interior, Tij its boundary consisting of linear segments Tij^i±ij ^rij^ij±i. 
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+ ^Ni-lj^Ni-lj 
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+ ^iN2-l^iN2-l 
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= 0,-- 
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for i = 1 , • • • A^i — 1 and j = 1 , • • • A^2 — 1 , and 



^iN2j^iN2 

The elements A^j are defined in Appendix (see ([I9j)-([20j) and ([2l])^([22j) ). 

At each time level kr we can represent a solution u^j by a long vector 
u = by mapping the node (i, j) of the two dimensional spatial domain to 
the one-dimensional vector, i.e. / = = j - Ni + i for i = 0, . . . , A^i and 

j = 0, . . . , A^2 and setting Uj = u^j. System of linear equations ([6j) can be 
then rewritten as a linear equation 

Au = b (7) 

for the solution vector u = u^. The dimension of the square matrix A as 
well as of the vectors u and b is (A^i + 1)(A^2 + !)• 
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We can solve the aforementioned linear problem ([7j) by means of the 
Successive over-relaxation method (SOR) In this linear equation solver we 
iteratively solve 



(8) 



for / = 1, • • • , (A^i + 1)(A^2 + !)• We repeat the iterative process for p = 
0, 1, • • • ^Pmax until the prescribed tolerance level for the difference ||u^^+-^^ — 
u^^^ll < tol is achieved. The iterative SOR scheme (|8j) can be written in the 
operator form: 

^X^u(^) + 6^, (9) 

where = (D+cjL)-1(D-cjU), = cj(D+cjL)-16. Here D, L, U stand 
for the diagonal, sub-diagonal and upper-diagonal parts of the matrix A, 
respectively. The parameter G (1, 2) is the so-called relaxation parameter. 
It can be used in order to speed up the convergence of u^^^ ^ u by minimizing 
the spectral norm of the operator T^. 



3. The constrained level-set method and its numerical approxima- 
tion 

3.1. Time- space continuous framework of the constrained level- set method 

In this section we introduce our new constrained level-set method for 
image segmentation. We suppose that there are two disjoint subdomains 
and Vtout of the domain Vt such that Vti^ is a subset of the interior of the 
segmented object and Vtout lies outside the segmented object. Furthermore, 
we suppose that there are two prescribed functions v^w G C (Q) with the 
property such that w < v in Q and < in and > in \ and 
w > in Qout and w < in Q\ Qout- Notice that any function u fulfilling 
w < u < V inQ must be negative in and positive in ^out- Its zero level-set 
contains the set in its interior and Qout in its exterior. 

Our purpose is to construct a solution u = u{x^t) such that it satisfies 
the level set equation ([2j) in the open region where w{x) < u{x^t) < v{x^t). 
Moreover, we require that w{x) < u{x^t) < v{x^t) for all x G G (0,T]. 
The reason why to prescribe range bounds on the level-set function u is to 
keep the set Vtin inside and Vtout outside the zero level set = {x, t^(x, t) = 0} 
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approaching the boundary of a segmented object when t ^ oo. To this end 
we consider the following partial differential inequality problem: 

Ut = foT {x,t) enx {0,T], s.t. w{x) <u{x,t) <v{x), 

Ut > Q\/ ' (^g^^^ foT {x,t) enx {0,T], s.t. u{x,t) =w{x), (10) 
Ut < QV • (^g^^^ ioT {x,t) enx {0,T], s.t. u{x,t) =v{x), 

d^u = at on, (11) 
u \t=o = Uini in Q. (12) 

It can be also viewed as the following variational inequality problem: find a 
solution ueJC C X where X = W^^\{0,T) : L\n)) n L\{0,T) : W^^\n)) 
such that 

{A{u)^ (j) — u) > 0^ for each G /C, 

where JC is the cone: /C = {^x G A', w{x) < u{x^t) < v{x)^ for a.e. x G 
n,t G (0,T)} (c.f. Brezis [4, Chapter 2]). The operator A : X ^ i-^((0,r) : 
H-\n)) is defined by 

Ut ^ f ^o^u 



and (., .) is the inner product in L^((0,T) : L^(f])), i.e. 

Jo Jn ^ V 
5.^. Time-space discretized framework of the constrained level-set method 

The discretization of ffT0l)"ffT2l) follows exactly from the discretization of 
the level set equation ([2])-(|4l) when taking into account the range bounds 
constraints w{x) < u{x^t) < v{x) for a.e. x G G (0,T). At each time 
step we have to construct a solution u to the following linear complementarity 
problem: 

(Au)/ = b/ for / s.t. wj < Uj < vj, 

(Au)/ > hi for / s.t. u/ = wi, (13) 

(Au)/ < b/ for / s.t. U/ = vj^ 
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for /= l,...,(iVi + l)(iV2 + l). 

In order to solve the linear complementarity problem f fT3l) we make use of 
the so-called Projected SOR method (PSOR) pjj,|7|] adopted for our problem. 

For each index / = 1, • • • , (A^i + 1)(A^2 + 1)5 we repeat the following up- 
dating of the vector u^^^ : 



uf^^^ = min{max{u?^+'\^,},^,}, (14) 

until the prescribed tolerance level for the difference Hu^^^-'^) — u^^^H < tol is 
achieved. Suppose that the matrix A has a positive diagonal, i.e. ajj > 
for each /. Assume that u^^^ ^ u as p ^ oc. Hence u(^+^) (1 — cj) u + 
(b - (L + U)u). It means 

u/ = min{max{[(l - cj) u + cjD"^ (b - (L + U)u)]/, wj}, vj}. (15) 
Clearly, 

for each index /. If the strict inequality Wj < Uj < Vj holds for some index 
/, then 

uj = [(1 - cj) u + uB-^ (b - (L + U)u)]/. 

It means that 

(Au), = b,. 

On the other hand, if uj = wj then uj < vj. According to (fT31) we conclude 
that 

u/ > [(1 - a;) u + ul)-^ (b - (L + U)u)]/. 
Since a// > we obtain 

(Au), > b,. 

Analogously, if U/ = for some index / then it follows from (ITSll that 
u/ < [(1 - c^) u + cjD"^ (b - (L + U)u)]7. 

Thus 

(Au), < b,. 

Hence the vector u is a solution to the linear complementarity problem fITBll . 
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4. Application to image segmentation and computational results 

In this section, we present experimental results obtained by the con- 
strained level-set method. We first demonstrate the eflFect of the obstacle 
constraints applied to artificial images. Fig. [2] a) (left) shows an image we 
want to segment. There are two rectangles with centers at points (0.4, 0.5) 
and (0.6,0.5). The width of the rectangles is 0.1 and the height is 0.4. The 
thickness of the rectangles edges is 0.04. On the inner edge of each rectangle, 
there is a thin hole. The image intensity function is defined on the domain 
= (0, 1)^. The initial curve is a circle with the radius r = VO.08 centered 
exactly between the rectangles. Its signed distance function (taken as the 
initial level-set function Uini) is depicted on Fig. [2] a) (right). The numerical 
mesh consisted of 128 x 128 grid points. If the thin holes in the rectangles 
should be taken into account we may want to segment only the rectangles 
edges. It can be achieved by setting the regularizing parameter 6 = 1. In |12| 
Mikula and Sarti interpreted the regularization parameter e as a parameter 
by which we can control "convexity" of the final segmentation curve. The 
result taken at time t = 1.2 is depicted at Fig. [2]b). On the other hand, if 
small holes are present by mistake, we may want the segmentation curve to 
fill them. It may be achieved by setting e = 0.0001. However, a convex hull 
of both rectangles is segmented as we can see on the Fig. [2]c). 

The segmentation we aim to can be achieved by prescribing one constraint 
guaranteeing that the part of the image between the rectangles must be 
outside the segmentation domain. We construct the constraint by putting a 
red bar placed on between the rectangles on the Fig. [3] a). Its width equals 
0.04 and the height is set to 0.6. The constraint function w{x) is positive on 
the red bar region and negative everywhere else. We again set e = 0.0001. 
The result taken at time t = 67 is shown on Fig. [3]b). 

In our last synthetic example shown in Fig. H] we place one more con- 
straint inside the left rectangle. The segmentation poses features of both 
experiments on the Figures [2] and [3l It shows that the constrained level-set 
method is capable of controlling which part of the image we want to segment. 

In the remaining two figures we present a real application of the con- 
strained level-set method to segmentation of cardiac MRI data. First we 
want to segment both, the left and the right ventricle of a heart. The usual 
level-set method for the image segmentation (c.f. [12]) is not capable to sep- 
arate them. In Figure [5] we present a) the initial set-up of the constraint 
barriers; b) segmentation without obstacles involved; c) segmentation with 
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c) 



Figure 2: a): We plot an initial curve set as a circle with the radius r = V0.08 (left) 
and its level-set function (right), b): The segmentation at time t = 1.2 after setting the 
regularizing parameter e from ([5j) to 1 is depicted together with the level-set function, c): 
the segmentation result with e = 0.0001 at time t = 145. 
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b) 

Figure 3: a): An example of prescription of a one constraint depicted as the red bar in the 
middle of the image. The initial curve is a circle with the radius r = VO.08. b): The 
segmentation at time t = 67 shows nice separation of both rectangles. The corresponding 
level-set functions are depicted on the right. 
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b) 

Figure 4: a): An example with two prescribed constraints depicted by red bars (top left). 
The initial curve is a circle with the radius r = V0.08. b): The segmentation result 
at time t = 67 is shown in the bottom left part of the figure. The corresponding level-set 
functions are plotted on the right. 
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a) b) c) 

Figure 5: a) Setup for the constrained level set method with the initial curve (green) 
and two constraints - the red curve stands for the exterior while the blue one is for interior 
of the segmented region, b) Segmentation results obtained by the unconstrained level-set 
method, c) Segmentation obtained by means of the constrained level-set method. The 
interior constraint helps to capture the bottom of the right ventricle. 

constraints separating both ventricles well. Figure [6] demonstrates more com- 
plex segmentation result in which we separated the left ventricle, septum and 
the right ventricle together with pericardial fat. 

Conclusion 

We proposed a new constrained level-set method for the image segmenta- 
tion. The method gives a possibility to an expert to prescribe an additional 
information concerning the expected shape of the segmented object. In this 
method, we may preset fixed constraints telling us which parts of the seg- 
mented object should be inevitably inside the segmented region and some 
parts must remain outside. The complementary finite-volume method was 
used for the numerical approximation of the level-set equation. It was com- 
bined with the Projected Successive Over Relaxation method for solving the 
corresponding parabolic variational inequality problem. We demonstrated 
the diflFerence between the new and the usual level-set method on several 
artificial images as well as on data from the magnetic resonance. 

Appendix 

In this section we show details of derivation of the full time-space dis- 
cretization of (Ej). As usual in the finite volume methods, we apply the 
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a) 



b) 



Figure 6: a) Setup for the constrained level-set method, b) Segmentation separating the 
left ventricle (small almost circular curve inside) surrounded by septum (larger curve) and 
the right ventricle with pericardial fat (the outer curve). 



Stokes theorem to the right hand side of (Ej). We obtain 

where v is the outer unit normal of Vij . We have to approximate the " fluxes" 
/r - g^^^"^ "capacity" numerically. We assume that the 

"fluxes" are constant on each segment F^^ - of the finite volume boundary 
Tij. It means that 



Qk-1 y "^3.^3 Qk-\ ^ "^i^id "^3.^3- 

First we evaluate the norm of smoothed gradient of the image intensity 
function Sij = |Gcr * V/o|^^-. Here Gcr{x) = (27r(7^)"-^ exp(— ||x|p/2(7^) is the 
Gaussian mollifier. We refer the reader to [12] for details how to effectively 

calculate the convolution G(j * V/q = VGcr * ^o- Then g^- = 1/ ^1 + Xs^- on 

ujh and we approximate g^. .-. on the finite volume edges as follows: 

Iqo Iqo 

9ij,i±ij = 2 ^^^3 ~^ 9i±ij) 5 9ij,ij±i = 2 ^^^3 ~^ 9ij±i) ' 
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Denote Vu^. -• = idx.u^- --.dxnU^- - ]. The approximation of \/u^. .-. in the 
direction of the vector u^j fj is obvious: 

In order to calculate remaining coordinates of Vu^j - which are perpendicular 
to u^j fj we need to know the value at the ends of /^^ - (corners of the finite 
volume Vij). They can be only approximated using the value of at the 
neighboring volumes: 

where p = i±l^q = j±l. Hence we obtain the approximation of Vu^j - in 
the direction perpendicular to u^j fj in the form 

(18) 

Having calculated approximation of Vu^.^. we can approximate the term 
Q^. .-. as follows: 



where ]9 = i±l,g = j''±l. The "capacity" term and the function 

are approximated by a constant value on the finite volume Vij. Taking the 
averaged value 

yields the approximation 

1 - , , 1 4- - ^ ^ 1 4- - 

It follows from ( IT6l) that 
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Notice that u^j fj = {i — i^j — j) (see Fig. [T]). In such a regular grid, one 
coordinate of u is always vanishing. It cancels one coordinate of Vu^ in the 
inner product Vu^ • Recall that /^^ - attains only the values hi or /12. We 
have 

,k \ ri^ n/k n.k 



Vzi \ 9ij,i+ij "ij , , 9ij,ij+i "ij+i u 



-\-tl2 r h ^173—1 T 



Finally, equation ( ITOl) is approximated as 




ij,ij±l 



Let us denote 

^0 ^0 

(19) 
(20) 

To approximate the Neumann boundary condition Vu - u = on dQ where 
u is the outer unit normal of the domain Q we set 



hi hi 

- 4 _ 0, ^-^^ ~ ^-^^-^ = for z = 0, . . . , iVi , 



h2 h2 
which yields 

= 1, A'l^ = -1, A^^^. = 1 and = -1 for J = 0, • ■ • , N^, (21) 

A^o = 1, = -1, A^^^ = 1 and = -1 for i = 0, • ■ • , iVi. (22) 
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